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Abstract 


In this paper we introduce a new class of stochastic Petri nets in which one or more places can hold fluid rather 
than discrete tokens. We define a class of fluid stochastic Petri nets in such a way that the discrete and continuous 
portions may affect each other. Following this definition we provide equations for their transient and steady-state 
behavior. We present several examples showing the utility of the construct in communication network modeling 
and reliability analysis, and discuss important special cases. We then discuss numerical methods for computing the 
transient behavior of such nets. Finally, some numerical examples are presented. 
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1 Introduction 


One of the difficulties encountered while using Petri nets is that the reachability graph tends to be very large 
in practical problems. Drawing a parallel with fluid flow approximations in performance analysis of queueing 
systems, we may define fluid within a Petri-net to approximate token movement. Alternatively, some physical 
systems have not previously admitted a Petri net modeling approach, as they explicitly contain continuous 
fluid-like quantities which are controlled with discrete logic. This paper presents a new methodology for 
modeling such systems. 

Stochastic fluid flow models are increasingly used in the performance analysis of communications [3, 10, 13] 
and manufacturing systems. On the other hand, stochastic Petri nets with discrete places provide a useful 
framework for specifying and solving performance and reliability models of discrete event dynamic systems 
[1, 6, 9, 17, 19]. It is natural to extend the stochastic Petri net framework to Fluid Stochastic Petri Nets 
(FSPNs) by introducing places with continuous tokens and arcs with fluid flow so as to handle stochastic 
fluid flow systems. This paper extends the model in an earlier paper [18] by allowing the level of fluid 
in continuous places to affect the enabling of timed transitions and the rates of fluid flow into and out 
of continuous places. Rules for transition enabling and firing are extended to reflect the notion that flow 
through a fluid place represents token movement. 

An FSPN contains two types of places: discrete places containing a non-negative integer number of 
tokens, and continuous places containing fluid. Transition firings are determined by both discrete places 
and continuous places, and fluid flow is permitted through the enabled timed transitions in the Petri net. 
Associating exponentially distributed or zero firing time with transitions, we can then write the differential 
equations for the underlying stochastic process. We also provide additional examples of FSPN usage, and 
discuss numerical issues arising in the solution of the underlying dynamic equations. 

The main motivation of this paper is to put the research by Mitra and his colleagues [3, 10, 13] in the 
context of Petri nets, to make some extensions and to investigate the numerical transient analysis of such 
stochastic fluid models. 

The paper is organized as follows. In Section 2 we develop the fluid model of stochastic Petri nets and in 
Section 3 we discuss their analysis. Examples and special cases are described in Section 4, numerical solution 
techniques are described in Section 5, while numerical examples are given in Section 6. The paper concludes 
in Section 7. 

2 The Stochastic Fluid Model 

Following the customary notation [8, 14, 16] for defining Petri nets and their extensions, we define a fluid 
stochastic Petri net (FSPN) as a 8-tuple (V, T, A, m 0j T, W, ft, Q ). V is a set of places partitioned into a 
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set of discrete places Vd and a set of continuous places V c . The number of fluid places is F > 0, indexed by 
k = 1, 2, • • • , F. In a graphical representation, we shall depict continuous places by means of two concentric 
circles. The set of transitions T is partitioned into a set of (exponentially distributed firing) timed transitions 
Te and a set of immediate transitions 7 } . The set of directed arcs A is partitioned into two subsets A c and 
Ad- A c is a subset of (V c x Te) U (Te x V c ) while Ad is a subset of (Vd xT)U(T x Vd )• In a graphical 
representation, arcs in A c are drawn as double lines (to suggest a pipe) while those in Ad are drawn as single 
lines. 

Let rrid = * € Vd) be the vector of the number of tokens in discrete places and let x — (x*, k £ V c ) 

be the vector of the fluid levels in continuous places. We will say that rrid is the discrete marking of the 
net. Let M denote the number of discrete markings, which will be indexed by the symbols i and j, with 
i,j = 1, 2, ♦ • M . The complete state (marking) of a fluid Petri net is described by the vector m = (x, m<f) 
where is the marking of the discrete part of the state and x keeps track of the fluid levels in the continuous 
places. Let M be the set of all complete markings (x, m^) and Aid be the set of all discrete markings. The 
initial marking is mo = (xo, mao)- 

In our formulation an enabled transition in Te may drain fluid out of its continuous input places, and 
may pump fluid into its continuous output places. The rates of flow may be dependent on the complete 
marking (x, m^). In a general formulation of a stochastic Petri net (embodied, for example, by SPNP [7, 8]) 
the conditions for enabling can be specified either through explicit arcs or through Boolean functions known 
as guards. We will allow both of these possibilities. We will continue to use the enabling and firing rules 
employed in SPNP with the additional possibility of a guard associated with a timed transition being able 
to base the enabling condition not only on the discrete marking of the net but also on the continuous part. 
We disallow the enabling of immediate transitions by fluid levels as this would lead to x-dependent vanishing 
markings, which cannot be eliminated in a manner analogous to that of GSPNs. Thus the guard function Q 
is defined for any timed transition in Te so that Q : T x M — ► {0, 1}. For a timed transition t £Te, Q(t , m) 
is a Boolean function that will be evaluated in each marking, and if it evaluates to true , the transition 
may be enabled; otherwise r is disabled. Upon firing, the transition removes a specified number of tokens 
from each discrete input place, and deposits a specified number of tokens in each discrete output place. 
The basic extension we have made from [18] here is that fluid levels in continuous places can change the 
enabling/disabling of timed transitions in the discrete part of the net. A discrete marking rrid is said to be a 
vanishing marking if one or more immediate transitions are enabled by it; otherwise it is a tangible marking. 

The firing rate function T is defined for timed transitions Te so that T : T e x M — ► 1R + . Thus if a 
timed transition r is enabled in (tangible) marking m, it fires with rate T(t, m). Note once again that in 
[18], these rates were not allowed to depend upon the fluid levels but now we do allow the firing rates to be 
dependent on fluid levels. 

As in [18], the weight function W is defined for immediate transitions Tj so that W : T j x Md — ► IR + . 
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Thus if an immediate transition i is enabled in (a vanishing) marking m d , it fires with probability 

W(i, m d ) 

T: W(0,m d ) 

$£Ti enabled in m& 

Next we describe the evolution of the continuous part of the marking. The flow rate function 7 Z is defined 
for the arcs connecting a continuous place and a timed transition so that H : A c x M — ► IR + U {0}. Thus 
when the FSPN marking is m t € M at time t, fluid can leave place k G V c along the arc ( k , r) G A c at rate 
H((k, r), m t ) and can enter the continuous place k at rate 7£((r, k), m t ) along the arc (r, k) G A c for each 
T g Te that is enabled in m t . The instantaneous rate at which fluid builds in a place k G V c at time t, in 
marking m tl is then given by 

E r), m t ). 

t£T e enabled in m 1 t£Te enabled in m ( 

We require that for every discrete marking m d and arc (r, k), the rate k), (f, m d )) be a “nice’’ function 

of x, e.g.. it is piecewise continuous. Observe that since m t contains continuous levels x, rk{m t ) may change 
as a function oft even if the discrete part of m t does not change. Once again we have extended the definition 
in [18] by allowing these rates to be dependent on the fluid levels. 

Now let A'jfc(t) be the fluid level at time t in a continuous place k G "P c - We assume that there is an upper 
bound on the fluid content, that is, X k (<) < B k for all t > 0. If there is no such upper bound, we set B k to 
oo. Then the sample path of X k (t) satisfies the differential equation 

' [r*(m ( )] + if X k (t) = 0 

dX k (t) _ if X k (t) = B k 

dt ~ r k (m t ) if 0 < X k (t) < B k and r fc (m < _)rt(m <+ ) > 0 

0 if 0 < X k {t) < B k and r t (m < _)r i (m t+ ) < 0 

In the case X k {t) = 0 and r k (m t ) < 0, we set the actual rate equal to zero (denoted by [nt(m t )] + = 
max(rjt(m t ),0)) in order to maintain X k (t) > 0. In the case that X k {t) = B k and r k (m t ) > 0, we set the 
actual rate equal to zero (denoted by [r*(m t )] _ = min( r k (m t ), 0)) in order to maintain X k (t) < B k . For the 
explantion of the remaining cases, we refer the reader to [10], Section II. The key observation (for the fourth 
case) is that a sign change from + to - in r k (m t ) at m, will “trap” X k (t) to be constant. Finally, let M d (t) 
be the discrete marking at time t. 

In the next section we study the joint process (A'(f), M d (t)), where A"(t) = [X k (t),k G V c ]- 
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3 Analysis 


Recall that A'*(<) is the fluid level in the k th continuous place at time t. The reachability graph corresponding 
to the discrete part of the net gives rise to a stochastic process that is a Markov process with the state space 
M. In the special case that the discrete part of the net is not affected by the fluid levels, discrete part of 
the net gives rise to a CTMC[1]. Let S be the discrete state space and let Q(x) = [<Hj(%)] be the matrix 
of transition rates derived from the firing rate function T of section 2. S corresponds to the set of tangible 
discrete markings in M d • For every x and place k £ V c define the diagonal matrix Rjt(x) = diag(rk(x , m^)), 
m d E S ). 

Define the distribution function H(t,x,m d ) = P{X(t)<x , and M d (t) = m d ) and let H{t,x) = 
[. H(t , f, rrid ), rrid E 5] be a row vector. 

To begin with, we assume that the rate functions Rk{x) are differentiable functions of x and transition 
rates Q(f) are piecewise right continuous functions of x . We also assume that the capacities of the continuous 
places are infinite. Under these assumptions, it can be shown that (X(t), M d (t)) has a density h(t,x,m d ) 
for all x with non-negative components and m d £ S and has a probability mass c(t,x, m^) if x has at least 
one component equal to 0, where 

u - x _ P(Xk(t) = 0 if Xk = 0, X k (t) E (x k ,x k + dx k ) if x k > 0 Vi) 

C\t j X , 77l d ) — -1—r , 

1 iJt Xlc >0 

Let h(t,x) denote the corresponding row vector of m^). 

The next theorem gives the coupled system of partial differential equations satisfied by h and c. These 
equations describe the transient behavior of the FSPN. 

Theorem 1: 

The equations satified by h are 


dh 5(/iRfc(£)) 

dt “ Q Xk 

k 

& 

T-C 

II 

(2) 

c(t,x, m d ) 

= 0 


—c{t,x,m d ) + h{t,x,m d )r k (x,m d )+ 

k:x k = 0 

if, for any k, x k = 0 and r k (x, m d ) > 0 

(3) 

k: i k >0 ° Xk 

= ^2c{t,x,i)q itmi (x) 
<es 



if Xk = 0 => r*(£, m<i) < 0 Vi (4) 
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Proof: 


A rigorous proof can be obtained by using the same technique as in [12]. Here we provide an intuitive 
proof for the ^-equation. 

Assume for simplicity only one fluid place and ri(x,i) > 0. Consider a portion of the function h(t,x,i) 
at one instant in time in the vicinity of the location x = x/. We consider a cell surrounding this point whose 
left and right boundaries are located at x_ and x+ respectively. Let the discrete values h(t, xi,i) and qij(xi) 
represent the mean values of h and q t j respectively within the cell. This situation is illustrated in Figure 1. 



Ax 


Figure 1: Derivation of FSPN equations 

Probability mass must be conserved. We may therefore derive a balance equation for the change in 
probability mass inside the cell: 

change in probability total mass total mass 

mass in cell entering cell leaving cell 

Probability mass enters the cell at the location x_ at the rate r(x_, i) and leaves at x+ at the rate r(x + , i). 
In addition the cell gains probability from the corresponding cell of each other discrete marking m j , j ^ i 
at the rate qji{x\) and loses probability to them at a total rate of — qu{x\ ). The balance equation for the cell 
probability mass becomes 

Az d/t (*’ X ~ ) = h(t,x-,i)ri(x-,i) - h(t,x + ,i)r 1 (x + ,i) + Ax(^2q j i(xi)h(t,xi,j) + qi i (x,)h(t,x l ,i)). (5) 

j** 

The equation is identical for the case 7 * 1 ( 2 :, i) < 0. Dividing the vector form of equation (5) by Ax and taking 
the limit as Ax — * 0 we obtain 

dh d{hRi(x)) _ r 
di + Tx " ftQ(x) • 

The argument generalizes easily to the case of more than one fluid place, yielding equation (2). The cell is 
in general a hypercube of sidelength Ax. 
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The boundary conditions follow from the fact that h(t,x, i ) is a probability density function. 


□ 


Theorem 2: 

The equations for the cumulative probability distributions H are 

# + E &*••(*» - (6) 

with the boundary conditions 

H(t,x,i) = 0 at Xk = 0, if r*(x, i) > 0 (7) 

3H 

lim — = 0. (8) 

Tk — oo OXk 


Proof: 

Using the abbreviations 


d JL 

dx 


d^H 

dx\ . . . dxp 


[ dx = [ F . . . [ 

Jo Jo Jo 


dx i . . . dxp 


and noting that 


H(t, x, i) = c(t, £, i) + / 

Jo 

we substitute into Equation (2) and integrate with respect to x\ 


l 


d dH \ „ 
-~ ) dx + 

0 ut uX 


jf? 


5 dH n 

a5T a?®"' 1 ' 


) df = i 


which yields (6). Equation (4) is implicitly contained in Equation (6), along with the boundary condition 

( 7 )- 


The boundary condition (7) follows from the observation that the fluid level in place k cannot remain 
at 0 for a positive amount of time in state m<f if the net rate r*(x, m<*) > 0. The boundary condition (8) 
comes from the observation that H represents a probability distribution with respect to each xjt, and must 
therefore approach an asymptotic value Xk tends to infinity (or reaches its maximal value). 


□ 

This method of deriving partial differential equations that represent conservation laws is well established 
in computational fluid dynamics, where it is known as the “Finite Volume” technique [15]. 
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The assumptions on R*(x) can be relaxed and we can allow R*(x) to be piecewise differentiable and 
fluid places to have finite capacities. This may introduce non-zero probability masses at the inter-region 
boundaries and will need to be explicitly accounted for. 

The domain of Equations (2) and (6) is 


0 < x k < B k 
0 < t < oo, 

although in practice we will only be interested in the finite domain 0 < t < tmax , 0 < x* < min {xmax, Bk) 
(where xmax is finite when Bk — oo). 

The initial conditions for Equation (6) are 

H( 0, x, i) = 1 if * = mdo and x > x 0 

H(0,x, i) = 0 otherwise. 

The initial conditions for Equation (2) are 

c(0, x, i) = $(m 0 ) 


where 6 is the delta function. 


Now suppose the following limits exist 


F(x) — lim H(t, x). 

t — ►00 


Then from Theorem 2. we see that the steady-state distribution F(x) obeys the following system of differential 
equations 

s h, (^> r ‘ w ) = - 1; ^ m 

with the normalization condition lim F(x)e = 1, where e is a column vector of all l’s. 


We note that the steady-state distribution F exists when 


lim V ?Tt(x)rjfc(x, i) < 0 = 

t 


where 7r(x) is the solution of 


7t(x)Q(x) = 0 

= 1 
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3.1 FSPN with a Single Continuous Place 

In the special case of a single continuous place, Equation (9) reduces to: 

±(F(x)R(x)) = F(x) Q(ar) - jT dt. (10) 

In the following subsections, we consider three special cases of the Equation (10). 


3.1.1 Constant Case 

In the special case that R(x) and Q(x) are both independent of x, following [3], solution of such an equation 
is of the form F(x) = he Xx where A is a row vector and A is a scalar. Substituting in (10) we have 


A(AR-Q) = 0. (11) 

If a non-zero h is to satisfy the above equation, we must have det(AR — Q) = 0. The number of solutions 
of det(AR — Q) = 0 equals the number of non-zero diagonal elements of R(x). Let these solutions be denoted 
by Aj, A 2 * , Afc. Let A t be the solution to A,(AjR — Q) = 0. Then the general solution to (10) is given by 

k 

F(z) = ^a i h,e A ' r (12) 

i = l 

where the scalars a* need to be determined from the boundary conditions and the boundedness of F(x). It 
is known that the number of A;’s with positive real part equals the number of negative diagonal entries of 
R. The coefficient a, corresponding to an eigenvalue A* with Re(Xi) > 0 must be zero in order to maintain 
boundedness of F(x). The remaining coefficients a,* are uniquely determined by the boundary condition 
F( 0, m) = 0 if r(m) > 0. 

Now if R(x) and Q(x) are both piecewise constant functions of x, we can apply the above procedure 
for each different segment and piece the individual solutions together [10]. In the general case, we can use 
numerical solution methods for linear odes that are available. We refer here to explicit methods such as 
RKF-45 or implicit methods such as implicit Runge-Kutta [4] or TR-BDF2 [5] and so on. 


4 Examples 

Next we examine a number of examples to illustrate the modeling power of FSPNs. 
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c(x) 



:> 


■> 


W P V #P 2 ,x) 


L(#p x ,#p 2 ,x)U 



Figure 2: FSPN Model of Statistically Multiplexed Network Switch 


4.1 Statistical Multiplexing of a Network Switch 

For our first example, we recast the problem studied in [10] as a FSPN. We have K sources of ATM cells. 
Each source emits both high priority and low priority cells; the arrival rate of each depends on the state of 
a two-state Markov chain, e.g., if the source’s Markov chain is in state s € [1, 2], then the rate at which high 
priority cells is produced is A^\ and the rate at which low priority cells is produced is A|^. All cells are 
delivered to a single switch, with buffer capacity B. All low priority cells are discarded when the buffer level 
is greater than some B t < B. 

The FSPN for this problem is illustrated in Figure 2. The number of tokens in place pi reflects the 
number of sources in state i , i = 1 ? 2. The rate at which sources change from state 1 to 2 (alt., 2 to 1) 
is #piAi (alt., #p 2 A 2 )• Letting r* tgh and rj ou/ denote the rate at which high and low priority cells are 
generated by a source in state i , the aggregate arrival rate of fluid to the buffer from sources in state i is a 
function of the discrete marking : 

/i(#Pl» #P 2 ) = #P>(r>' 9h + r' ou '). 

The overall arrival rate as a function of the discrete marking is 

T(#Pl,#P2) = /l(#Pl,#P2) + /2(#Pl,#P2)- 

The fluid level in the place loss reflects the total cell loss since time 0 (like the cell buffer, it is initially 
empty). The rate c(x) at which cells are successfully switched out of the buffer with level x is c(x) = c if 
x > 0, and c(x) = 0 if x = 0. Of more interest is the function L(#p\,#p 2 ,x) describing the rate of cell loss 
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*r 

=]^-H 


Figure 3: FSPN for the Machine Breakdown Model 


(both into and out of the transition so labeled): 


£(#Pi,#P2,z) 


[t(#Pi > #P2) - c] + if 0 < X < B t 

#Pir‘r +#P 2 rl? w +[#pir h 1 i9h + #p 2 r h 2 igb -c]+ if B t < x < B 
#Pir‘r + *P 2 r l r + {#Pir\' 9h + #P 2 r h 2 i9h - c]+- 
. [c - #Piri' gh + #P 2 f 2 ,i,ft ] + if x = B t 


When x < Bt , low priority cells are not discarded automatically, and any loss is the difference between the 
aggregate arrival rate and the service rate. Low priority cells are dropped automatically when B t < x < B, 
and further loss may be due to an inability of the server to keep up with the aggregate high priority arrival 
rate. 


4.2 Machine Breakdown Model 

In the previous example, fluid levels have no effect on the behavior of the discrete portion of the FSPN. The 
next example shows how fluid levels may affect the firing rate of discrete transitions. Consider a system with 
jV statistically identical and independent components, each with a failure rate that depends on the overall 
system load [11]; the repair rate is fi. Work arrives to the system at a constant rate r and is completed at 
rate d per functioning machine. Work here is considered to be a non-negative real quantity. We model this 
system as an FSPN shown in Figure 3. 

The model has two discrete places ( p u and pd ), one continuous place and three timed transitions (tj , t r 
and The number of tokens in p u models the number of functioning machines (with the initial value 

N) while the number of tokens in pd is the number of failed machines undergoing repair. The firing of the 
transition tj represents a machine failure; given positive failure rates Ai and A 2 and positive scaling factor 
a, the load-dependent firing rate of this transition with fluid level x is taken to be 

f{*Pu,x) = #Pu(*i + A 2 (1.0-c-°* )). 


10 





Off 



Figure 4: FSPN model of an alternating switch 


Repair of a machine is modeled by the firing of t r ; its firing rate is p times the number of tokens in pd . The 
transition t a iways is always enabled and it continuously pumps fluid at rate r into the continuous place. The 
rate of firing of this transition can be chosen to be any positive value. Whenever transition tf is enabled, it 
drains fluid from the continuous place at rate d times the number of tokens in place 

4.3 Alternating Switch 

In the next example we model a communcation switch with two arrival streams. One stream is bursty, with 
a high arrival rate ruigh when active. It is described as an alternating process, that is on for an exponential 
period of time with rate A on , and off for an exponential period of time with rate A <>//. The other stream is 
slow, but constant (n ou ,). The switch services workload from either stream at rate p. 

The switch is designed to allocate exponentially distributed time-slots to each stream, with rates A high 
and A low . The time slots alternate— the fast stream is given a slot, then the slow stream, and so on. The 
FSPN modeling this switch is illustrated in Figure 4. 
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5 Numerical Transient Analysis 


5.1 Discretization 


We choose to solve the equations for the probability distribution (6) rather than the density equations 
(2), since the latter would require the numerical treatment of a delta function used to describe the initial 
condition. We describe here the case R(x) = const, Q(£) = const. 


We perform a semi-discretization of (6) in the coordinate directions using a first-order upwind method. 
In the nomenclature of section 3, the upwind scheme approximates values at the boundaries of the cell by 
those at the neighbouring grid points in the upstream direction, relative to the flow direction r. We choose a 
finite domain 0 < x* < xmax and use G equidistant grid points. The grid spacing (mesh size) Ax is therefore 
equal to xmax/(G — 1). Note that in order that the boundary condition (8) be approximately satisfied, 
xmax must be chosen to be sufficiently large. We use the notation = rk(l\ Ax, . . IpAz, i), 

- Ax, Ax) and = H(t, l\Ax } . . . , IpAx, i) . The upwind discretization is 

given by 


- — H(t, ...,x k = Ik Ax , ...,») ss 
ox k 


(i*i U,..(0 - Hi , (*)) if nt Ik j > 0 

(^» ( l+l) fc . ..(0 - ^ < 0 


(13) 


We obtain from the semi-discretization the linear system of ordinary differential equations 

dH 


dt 


= HQ 


(14) 


where the unknown vector 

^ 1 ,0 j , ...,0/? : 1 Af ,G — li,...,G— I f ^ 

is obtained by a lexicographic ordering of the unknowns by fluid place and by the discrete marking. 

The matrix Q is given by 

Q = Q + W 

where Q represents the CTMC of the discrete part of the net, and W the discretization of the space 
derivatives multiplied by the flow rates. 


The matrix Q is given by 


Q = 


Qn ■ ■ Qim 
Qmi ■ Qmm 


12 



where 


Qij 


Qij — 




Note that the matrix Q is a CTMC matrix. The block-diagonal matrix W has the form 


W = — 


1 

Ax 


W x 


0 


0 

w M 


where Wi G JR° FxGF represents the discretization of the second term of equation (6) for the z-th discrete 
marking using the upwind method (13). Each Wi is a sparse, banded matrix in which each column is either 
zero or has a positive main diagonal coefficient and non-positive off-diagonal coefficients containing the rate 
values r { . In addition, the main diagonal coefficient is the negative sum of the other entries in its column. 
Each Wi, and thus also W, therefore represents the transpose of a CTMC. 

The linear system of ordinary differential equations (14) thus defined may be integrated by any standard 
method. In the numerical examples in section 6 we will use the Forward Euler scheme 

H{t + At) « H(t) + AtH{t)Q. 


Note that for this scheme the discretization mesh sizes must satisfy 

max rk At < Azr (15) 

in order that the integration be stable. 


5.2 Numerical Integration in a Transformed Domain 

In order to avoid the difficulty of solving equation (6) in an arbitrarily large and a priori unknown domain, 
we can perform the coordinate transformation 

Vk = 1 ~ (16) 

which maps the infinite interval x* = [0,oo] to the finite interval yt = [0, 1]. Equation (6) then becomes 

f + E &■*«<» - »)■ = - J’ < 17 > 

for The boundary and initial conditions remain unchanged. 

This minor modification to the equations makes their solution significantly easier, since the decision 
where to place xmax must no longer be made, and the danger of an inappropriate choice is avoided. The 
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coordinate transformation has the additional advantage of compressing the infinite interval of the large x 
values, where in many cases the solution shows virtually no structure, into a small space. Furthermore, 
we can obtain numerical values at x — oo, which correspond to the probabilities for the discrete markings, 
whereas in the untransformed case these must be approximated by values obtained at x = xmax. 

5.3 Problem Size 

Recall that G denotes the number of gridpoint.s in each dimension x x , F the number of fluid places, and M 
the number of discrete markings. The number of time-steps to be integrated is T . 

The computational complexity for the solution is 0(TMG F ) floating point operations, since for each of 
T timesteps we must increment each solution value in an F-dimensional grid of sidelength G for each of 
M discrete markings using 0(1) operations. Note that for an explicit integration method such as Forward 
Euler, because the condition (15) must be satisfied, an increase in G must ultimately be accompanied by a 
proportionate increase in T. 

The storage requirements of the algorithm are at least 4 MG F bytes since for each of M discrete markings 
we must store an F-dimensional grid of floating point numbers with sidelength G. Solutions at successive 
time-steps can be overwritten. 

Since the sidelength of the grids G may typically be of the order of 50 or more when a simple discretization 
is used, we see that this would seriously limit the size of the FSPNs that can be solved. Future work must 
therefore include strategies for reducing the amount of memory needed to represent the function H . 

5.4 Time Integration by Randomization 

Randomization is a numerical method widely used for the solution of systems of ODEs of the form (14). It 
has the advantages of high numerical stability and low roundoff error, a priori specification of absolute error 
tolerance requirements, and is in addition often found to be faster than numerical integration schemes. The 
method has superior roundoff error behavior when the matrix P = (J + ^Q) has all non-negative entries, 
where A > max|Q tt |. This is, for example, the case for a CTMC. 

As the following Theorem shows, the semi-discretized FSPN equation (14) also has this property, indi- 
cating that randomization may be the method of choice for computing transient solutions. 

Theorem 3: 

For the matrix P = (I + jQ) where Q is the matrix of equation (14) holds 

Pjj > 0 Vi, j 
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Figure 5: Numerical Solution of Breakdown Model 


Proof: 

It is sufficient to show 

Q<*,Wii<0, and Qy.Wy > 0 (18) 

“Q ii ^ Qj; i Qji ? ~f~ j (19) 

We have 

Q«,Wi< <0, Qtji W t j > 0 i?j 

q« = -X>, w« = -E w « 

»?*? 

since Q is a CTMC and W represents an upwind discretization. Equations (18) and (19) follow directly 
from Q = Q + W. 


6 Numerical Examples 

6.1 Machine Breakdown Model 

For our first illustration of the behaviour of an FSPN we choose the model of section 4.2. We consider the 
case of one processor only, and parameter values ofAi = 2, A 2 = 0, jx = 3, r = 1 and d = 2. In this case, 
both R and Q are independent of x. We solve the equations for the distribution (6) in the range 0 < t < 4, 
0 < x < 4. We discretized with stepsizes Ax = 1/64 and At = l.Oe - 4. The numerical results for this 
problem are shown in Figure 5. 

For this simple case, it can be shown via Laplace transforms that the steady state solution is given by 

H(x,UP) = 0.4(1 - e~ r ) 


15 



Figure 6: Numerical Solution of Modified Breakdown Model 


H(x, DOWN) = 0.2(3 - 2e" x ) 

The transient solution values obtained at t — 4 agree closely with these results. 

Now we modify the example to allow dependency of the matrices R and Q on the fluid level X. First we 
extend the model to contain a backup server which is only used when the primary server is down and the 
load level reaches a value of 0.5. We thus have 

{ 2 m d — UP 

0 m d = DOWN, x < 0.5 

2 m d = DOWN , x > 0.5 

The dependency of Q on x is obtained by setting 

A = \ x + \ 2 (l- e - ax ) 

choosing a = 1, Ai = 1 and A 2 = 2. The results of this computation are depicted in Figure 6. Note the 
discontinuity and change of slope at x = 0.5, when the backup server is started. 

Figure 7 shows the results for the unmodified breakdown model, using the equation in transformed 
coordinates (17). Here integration until t — 24 has been performed. Note that the solution obtained for 
H(DOWN) at t = 24 appears as a straight line through the origin with gradient 0.4. corresponding to the 
analytic solution of Equation (20). 

6.2 ATM Switch Model 

We now consider the multiplexed network switch model considered in section 4.1. We set the number of 
sources K to one, parameters B 1 = 3000 and B = 6000. High priority packets arrive at the rate 6 x 10 3 and 
low priority packets at a rate of 4 x 10 3 per second when the source is in state 1, and at rates 4 x 10 3 and 
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Figure 7: Probability distributions for breakdown model in transformed coordinates. Left: Server up; Right: 
Server down. 


Source state 1 




Figure 8: Probability distributions for ATM switch model. Left: Source state 1; Right: Source state 2. 


2 x 10 3 per second respectively when the source is in state 2. The switch can process both types of packet 
at a rate of 5 x 10 3 per second. The exponentially distributed rates of change between high and low priority 
packet generation are Ai = 0.5 and A 2 = 0.5. The results for the number of packets in the buffer are shown 
in Figure 8. The left and right results at t = 10 qualitatively match those of Elwalid and Mitra [10], Figure 
2, right and center, respectively. An appropriate choice of parameters also yields a result, not illustrated 
here, which is similar to [10], Figure 2, left. 


7 Conclusion 

We have defined a new class of stochastic Petri nets by introducing places with continuous tokens and arcs 
with fluid flows. This new class of fluid stochastic Petri nets (FSPNs) should be useful in modeling stochastic 
fluid flow systems, and may also be useful in modeling processes that control physical systems. Our model 
formulation permits the discrete and continuous parts to affect each other, endowing FSPNs with the ability 
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to both control the fluid flow, and have the discrete control decisions be affected by observed fluid flow. 
We have provided formal definition of FSPNs and developed the rules for their dynamic evolution. We 
have derived coupled systems of partial differential equations for the transient and the steady-state behavior 
of FSPNs. Spectral representation of the FSPNs with a single continuous place can be adapted from the 
literature on stochastic fluid flow models. We have presented a number of examples illustrating the modeling 
power of FSPNs, have considered issues arising in the numerical solution of the dynamical equations, and 
have provided numerically solved examples. 
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